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1.  Introduction 


Phononic  bandgap  materials — periodic  materials  that  block  acoustic  or  elastic  wave  propagation 
in  a  range  of  frequencies — have  a  variety  of  applications  including  piezoelectric  transducers  (1 ), 
ultrasound  and  polarization  filters,  acoustic  filters  for  sound  reduction,  and  vibration-less 
environments  for  sensitive  equipment  (2).  There  has  been  recent  research  into  designing  such 
materials,  first  dating  back  to  1994  (2),  which  showed  the  design  of  bandgap  materials  by  varying 
the  filling  fraction  of  a  specified  geometry.  Later,  the  actual  geometries  of  bandgap  material 
inclusions  were  optimized  (4).  Finally,  in  reference  5,  genetic  algorithms  were  used  to  maximize 
the  bandgap  of  phononic  materials  by  designing  an  inclusion. 

In  previous  work,  all  methods  used  a  discretized  unit  cell,  either  with  square  elements,  as  in 
reference  4,  or  triangular  elements,  as  in  reference  5.  A  major  downside  of  using  a  strict  grid  is 
that  the  grid  resolution  can  be  critical  in  outcome  and  efficiency.  Using  a  too  coarse  grid  limits 
the  complexity  of  the  solutions,  but  a  too  fine  grid  may  take  exceedingly  long  to  converge,  if  at 
all.  Another  issue  in  references  4  and  5  was  that  both  approaches  only  considered  two-phase 
systems.  While  not  a  severe  limitation,  it  may  be  fruitful  to  allow  the  optimizer  to  choose  which 
materials  to  use  in  a  design,  regardless  of  number.  Finally,  reference  5  only  considered  acoustic 
wave  propagation,  not  full  elastic  waves. 

These  limitations  are  addressed  here,  specifically  by  using  a  genetic  programming  approach  to 
geometry  optimization,  introduced  in  reference  6,  and  extended  in  references  7  and  8  and  using 
the  elastic  wave  equation  in  place  of  the  acoustic  equation.  The  genetic  programming  approach 
represents  geometry  using  a  tree  structure,  representing  combinations  of  arbitrary  convex 
polygons.  In  this  tree  structure,  each  leaf  node  contains  a  list  of  points  and  the  convex  hull  of  this 
list  of  points  represents  a  convex  polygon,  to  be  combined  with  other  polygons  as  defined  in  the 
function  nodes  of  the  tree.  The  lists  themselves  are  not  limited  in  size,  so  each  convex  polygon 
can  have  any  number  of  vertices.  This  approach  is  naturally  hierarchical  in  contrast  to  strict  grid 
methods,  as  in  references  4  and  5.  The  method  can  be  initialized  with  a  single  convex  polygon 
with  a  small  number  of  vertices,  while  more  polygons  and  vertices  can  be  inserted  later.  For 
multi-phase  systems,  material  parameters  can  be  inserted  in  the  leaf  nodes,  as  in  reference  8.  As 
is  shown,  lifting  these  limitations  leads  to  higher  quality  designs,  more  efficiently. 

The  remainder  of  this  report  describes  the  systematic  design  of  phononic  bandgap  materials  using 
genetic  programming.  Specifically,  section  2  gives  the  formulation  of  the  problem  in  terms  of 
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finite  elements  and  discusses  a  method  for  determining  the  bandgap  size  for  a  given  material. 
Section  3  describes  the  geometry  optimization  method,  which  uses  genetic  programming  and 
convex  polygons  to  construct  a  solution.  Several  optimization  results  are  presented  in  section  4 
and  the  conclusions  are  discussed  in  section  5. 


2.  Phononic  Bandgap  Material  Model 


The  main  purpose  of  this  work  is  to  lift  many  of  the  limitations  in  place  in  reference  5.  First, 
reference  5  only  considered  scalar  dilatational  waves;  here  the  full  vector,  elastic  problem  is 
treated.  Another  restriction  was  the  limit  to  two-phase  systems,  which  is  replaced  with  a  more 
general  method  capable  of  choosing  from  a  database  of  materials. 

Phononic  bandgap  materials  are  periodic  structures,  defined  by  some  unit  cell.  The  unit  cell  of 
the  problem  considered  is  shown  in  figure  1  and  has  a  few  defining  parameters:  the  side  lengths, 
li  =  1, | ,  and  the  acute  angle  a.  Figure  1  is  labeled  with  the  edge  vectors  1,.  though  typically  the 
problem  is  defined  using  side  lengths  and  acute  angle  via 


11  =  £ik, 

12  =  £2  cos(a)x  +  i 2  sin(a)y. 


(1) 


Figure  1.  Unit  cell  diagram. 


2 


In  a  phononic  bandgap  material,  assuming  it  is  infinitely  periodic  and  there  are  no  body  forces, 
the  displacement  u  =  uxx  +  uy y  follows  the  elastic  wave  equation 
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where  A  and  p  are  the  Lame  parameters.  Since  we  are  interested  in  time-harmonic  solutions, 
equation  2  can  be  converted  to  the  frequency  domain  assuming  ey(Jjt  time  dependence: 
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The  appropriate  boundary  conditions  must  now  be  formulated  to  reduce  the  problem  to  a  single 
unit  cell.  As  the  problem  is  infinitely  periodic,  we  assume  that  the  solution  will  not  grow  or 
shrink  as  it  travels  through  the  material.  (Were  it  to  grow  or  shrink,  it  would  either  grow  to 
infinity  or  shrink  to  zero,  neither  of  which  is  useful.)  Assuming  the  magnitude  does  not  change 
across  unit  cells,  the  solutions  on  opposing  boundaries  must  then  be  phase-shifted  versions,  also 
known  as  Floquet  boundary  conditions  (9),  which  can  be  written  as 

u(r  +  li)  =  u(r)ejk'h,  r  =  al2, 
u(r  +  12)  =  u(r)ejk''2,  r  =  61^ 

(4) 

V  •  u(r  +  li)  =  V  •  u(r)ejk'11,  r  =  al2, 

V  •  u(r  +  12)  =  V  •  u(r)ejk''2,  r  =  61i, 

where  0  <  a.b  <  1  arc  arbitrary  real  numbers. 

Equations  3  and  4  define  a  generalized  eigenvalue  problem  in  u  and  u,  given  a  wave- vector  k. 
Because  of  the  periodicity  of  the  problem,  solutions  will  occur  in  discrete  modes,  so  that  for  a 
given  k,  a  discrete  set  of  allowed  frequencies  of  propagation  u  can  be  found.  If  these  equations 
can  be  solved  for  all  k,  the  full  set  of  allowed  propagation  frequencies  can  be  recovered.  Regions 
of  no  allowed  propagation  can  exist,  which  defines  a  bandgap.  Fortunately,  the  solutions 
themselves  are  periodic  in  k,  so  only  a  finite  region  must  be  searched  for  bandgaps.  The  defining 
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region  of  k,  or  Brillouin  zone,  is  related  to  the  edge  vectors  1  by 


k,„  *  1m  —  2wSn 


(5) 


where  5mn  is  the  Kronecker  delta  function.  This  set  of  equations  can  now  be  discretized  using 
finite  elements  and  solved  as  an  eigenvalue  problem  to  find  the  size  of  the  bandgap  for  a  given 
phononic  bandgap  material.  This  process  is  described  in  the  following  subsections. 

2.1  Finite  Element  Formulation 


First,  the  governing  equations  are  converted  to  a  weak  form,  by  projecting  the  equations  onto  a  set 
of  testing  functions 
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where  txm  and  tym  are  the  x-  and  ^-components  of  a  vector  testing  function  tm  and  Q  denotes  the 
unit  cell.  The  double  derivatives  on  the  displacement  in  equation  6  can  be  removed  by  shifting 
the  derivatives  onto  the  testing  functions.  After  some  simplification,  we  get 
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Next,  the  unit  cell  is  meshed  with  triangular  elements  using  the  method  of  Persson  and  Strang 
(10).  One  important  note  regarding  the  meshing  is  that  the  nodes  on  opposing  boundary  edge 
pairs  should  be  equally  spaced  and  equal  in  number  to  facilitate  the  application  of  the  boundary 
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conditions.  Given  this  discretization,  the  solution  is  approximated  with  a  set  of  basis  functions, 
which  are  identical  to  the  testing  functions.  The  basis  and  testing  functions  used  here  are  linear, 
nodal  elements,  which  are  unity  at  a  node  and  zero  at  all  surrounding  mesh  edges.  (This  function 
is  also  known  as  a  hat  function.)  Expanding  the  solution  in  a  series  and  separating  the  boundary 
terms  gives 


N  M  P 

Mx(r)  «  J2vnbn( r)  +  ^2KK+n( r)  + 

n=  1  n=  1  n=  1 

N  M  P 

Uy(r)  «  J2Vnbn( r)  +J2Wnbn+N(r)  +  Znbh+N+M(r)i 
n= 1  n=l  n=l 


(8) 


where  ux  and  uy  are  the  x-  and  //-components  of  the  displacement,  &x  and  byn  are  the  x-  and 
//-components  of  a  vector  basis  function  bn,  N  is  the  number  of  internal  nodes,  M  is  the  number 
of  independent  boundary  nodes,  and  P  is  the  number  of  dependent  boundary  nodes.  The 
difference  between  independent  and  dependent  nodes  can  be  gleaned  from  equation  4:  A 
dependent  boundary  node  is  a  node  that  can  be  written  in  terms  of  the  solution  at  another 
boundary  node.  We  define  the  independent  boundary  nodes  to  be  those  nodes  along  li  and  12, 
except  at  the  r  =  li  and  r  =  12  comers.  Applying  the  Floquet  boundary  conditions  to  the 
unknown  boundary  term  coefficients  gives 


zp  =  eJ'kll<,  rp  =  rn  +  1,, 
zp  =  ejk'hw x,  rp  =  rn  +  12, 
zyp  =  ejk'hwyn,  r p  =  rn  +  1,, 
zy  =  ejk'hwyn ,  rp  =  rn  +  12, 


for  nodes  located  on  the  appropriate  boundary  according  to  equation  4.  These  relations  can  be 
collected  into  matrix  form 


zx  =  Dxwx, 
zy  =  Dywy. 


(10) 


The  matrices  Dx  and  Dy  contain  the  phase  shift  terms  from  equation  9,  one  coefficient  per 
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column.  After  substituting  the  approximate  solution  from  equation  8  into  equation  7,  we  get 
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Each  submatrix  in  the  matrices  above  correspond  to  specific  terms  of  equation  7  after  basis 
function  substitution.  The  terms  of  Bxx  equate  to 


Rxx  =  /  otx  bx  dr 

vv,mn  /  r°mwnK^  ? 


and  A**  is  given  by 
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The  dependent  boundary  terms  can  be  eliminated  by  substituting  equation  10  into  equation  11, 
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arriving  at 
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Equation  14  is  a  generalized  eigenvalue  problem  of  the  form 


orBu  =  Au. 


(15) 


As  in  reference  5,  ARPACK  is  used  to  solve  the  eigenvalue  problem.  Using  the  fact  that  B  is 
Hermitian,  positive  definite,  B  can  be  decomposed  using  a  Cholesky  factorization,  and  equation 
15  can  be  solved  more  efficiently.  Ultimately,  for  a  given  k  vector  (which  determines  the  matrix 
D),  the  eigenvalues  cu2  represent  the  allowed  frequencies  of  propagation  for  that  given  wave 
vector.  Knowledge  of  u  given  k  is  necessary  in  computing  the  bandgap  of  a  material  as  is 
discussed  in  the  next  subsection. 

2.2  Bandgap  Determination 

The  finite  element  formulation  is  only  the  first  step  of  determining  the  bandgap  of  a  structure  as 
any  single  solution  of  the  above  problem  only  reveals  the  propagation  frequencies  and  modes  at  a 
single  k  vector.  A  band  diagram  can  be  constructed  by  solving  the  problem  at  different  values  of 
k,  and  the  bandgap  can  be  approximated  visually.  Automating  this  process  requires  a  local 
search  method  to  determine  the  maximum  frequency  of  a  one  band  and  the  minimum  frequency 
of  the  next  highest  band.  In  reference  5,  the  bandgap  always  appeared  between  the  first  and 
second  bands  of  the  structure;  however,  because  they  are  composed  of  both  longitudinal  and  shear 
components,  elastic  waves  generate  bandgaps  between  higher  bands,  typically  the  third  and 
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fourth,  but  never  the  first  and  second.  Figure  2  shows  an  example  of  a  band  diagram:  Note  that 
the  bandgap  appears  between  the  third  and  fourth  bands.  Additionally,  the  band  structure  shown 
in  figure  2  has  mirror  symmetry  about  the  kx  =  ky  line. 


Figure  2.  Example  of  a  band  diagram,  generated  from  a  silicone  rubber/lead 
system  with  a  square  inclusion  in  a  square  unit  cell. 


Overall,  the  approach  is  similar  to  that  of  reference  5,  with  a  few  changes.  First,  the  relative 
bandgap  is  computed  by  finding  the  maximum  frequency  (versus  k)  in  the  first  band,  uq,  and  the 
minimum  frequency  in  the  second  band,  w2-  Given  these  quantities,  the  relative  bandgap  is 

7,  _W2-CUi 

"relative  , -  Oof 

\/W^2 


Previously,  a  simplex  method  was  used  to  determine  the  relative  bandgap.  The  simplex  method  is 
robust  and  fairly  efficient,  though  it  is  not  always  the  most  accurate  method.  (It  also  may  not  be 
more  efficient  than  a  derivative-based  method  with  a  good  starting  guess.)  Here  we  use  a 
standard  derivative-based  method,  the  Broyden-Fletcher-Goldfarb-Shanno  method  (BFGS)  (if). 
One  issue  with  derivative-based  methods  is  that  they  require  a  good  starting  location  to  find  the 
global  minimum  of  a  function.  As  can  be  seen  from  figure  2,  there  are  many  trouble  areas  for 
local  search  methods.  Instead  of  using  a  single,  random  starting  vector,  a  uniform  5-by-5  grid  of 
points  is  mapped  in  k-space  and  the  bandgap  is  evaluated  at  only  the  points  that  (inclusively)  fall 
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to  the  left  of  the  kx  =  ky  line.  The  best  initial  guess  of  the  15  points  is  used  as  the  starting  point 
for  BFGS,  and  initial  guesses  for  both  bands  of  interest  can  be  computed  simultaneously. 
Although  this  approach  requires  15  function  evaluations  initially,  it  typically  leads  to  fast  and 
robust  convergence  of  BFGS. 


3.  Geometry  Optimization  Method 


Phononic  bandgap  materials  may  be  composed  of  simply  a  matrix  with  a  single  inclusion,  or 
possibly  multiple  materials  each  with  arbitrary  topology.  The  method  outlined  in  reference  5 
considered  only  two-phase  systems  and  used  a  fixed  triangular  grid  to  represent  geometry.  The 
main  issue  with  the  approach  of  reference  5  was  the  fixed  discretization,  which  requires  more 
optimization  parameters  to  return  a  higher  resolution  solution.  An  increase  in  optimization 
parameters  leads  to  longer  run  times,  and  possibly,  poor  convergence.  Here,  a  more  flexible 
method  based  on  genetic  programming  first  discussed  in  reference  6  is  used  with  minor 
modifications.  The  main  advantage  over  a  fixed  discretization  is  that  large,  homogeneous  regions 
can  be  represented  with  far  fewer  parameters.  The  genetic  programming  approach  is  also 
naturally  hierarchical  because  points  and  polygons  can  be  added  to  chromosomes  as  needed. 

This  method  has  the  advantages  of  a  finer  discretization,  but  without  the  need  to  overestimate  the 
required  grid  size.  Subsection  3.1  reviews  the  geometry  encoding  and  subsection  3.2  discusses 
the  modification  to  the  genetic  operators  laid  out  in  references  6  and  7. 

3.1  Geometry  Encoding 

In  contrast  to  many  previous  geometry  encoding  schemes,  the  method  used  here  is  based  on  a 
genetic  programming  (72)  implementation  of  constructive  solid  geometry.  In  addition,  the 
primitives  used  to  generate  geometries  are  arbitrary  convex  polygons  as  defined  by  the  convex 
hull  of  a  list  of  points  (<5).  Some  methods  in  use  today  have  adopted  a  genetic  programming 
approach,  but  use  a  database  of  primitives  (13,  14).  In  subsection  3.1.1,  the  convex  polygon 
primitives  are  discussed,  and  in  subsection  3.1.2,  the  combination  scheme  that  accounts  for 
inhomogeneous  materials  is  reviewed  from  reference  8. 

3.1.1  Convex  Polygon  Primitives 

A  convex  polygon  can  be  defined  many  ways,  most  simply  as  a  polygon  in  which  any  line 
segment  with  both  end  points  interior  to  the  polygon  never  cross  a  boundary.  A  useful  way  to 
specify  a  convex  polygon  (for  our  purposes,  especially)  is  as  the  convex  hull  of  a  set  of  points. 
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The  convex  hull  of  a  set  of  points  can  be  defined  as  the  intersection  of  all  half-planes  that  contain 
those  points,  or  equivalently,  as 


Hn 


cPO 


k 

Xj  G  X,  >  0,  oti  =  l,  i 

i=  1 


(17) 


where  X  is  some  set  of  points  and  k  varies  from  1  to  the  total  number  of  points  in  X.  Several 
algorithms  exist  to  compute  the  convex  hull  of  a  list  of  points,  typically  these  algorithms  return 
the  subset  of  vertices  that  define  the  convex  hull.  There  are  different  algorithms  of  varying 
computational  complexity,  fortunately  0(n\ogn)  and  O (n  log  h)  (where  h  is  the  number  of 
vertices  on  the  convex  hull)  algorithms  exist  (15). 


Convex  polygons  are  ideal  in  some  ways  for  optimization.  First,  they  are  easy  to  define  using  an 
arbitrary  list  of  points  as  described  previously.  Second,  they  have  few  constraints,  i.e.,  any  list  of 
three  or  more,  non-collinear  points  (in  two  dimensions)  defines  a  convex  polygon  via  its  convex 
hull.  Third,  in  this  implementation,  the  number  of  points  in  a  list  is  not  restricted,  so  that  they  are 
flexible  in  “smoothness.”  Convex  polygons  can  range  from  simple,  triangular  shapes  or 
arbitrarily  “smooth”  approximations  of  truly  smooth  shapes.  (Of  course,  a  linear  interpolation 
scheme  is  used  throughout,  but  because  the  number  of  points  is  not  restricted,  they  can  approach 
smooth  shapes  such  as  circles  and  other  curves.)  Finally,  they  are  easy  to  manipulate  via  genetic 
operators,  especially  crossover,  as  is  discussed  below.  Using  a  database  method  can  be  difficult 
in  terms  of  crossover,  because  it  may  not  make  sense  to  combine  the  defining  parameters  of  a 
rectangle  (width  and  height)  with  a  circle  (radius). 

3.1.2  Inhomogeneous  Combinations 


While  convex  polygons  are  useful  for  optimization,  on  their  own  they  are  not  sufficient  to 
represent  any  arbitrary  geometry  or  topology.  Originally,  Boolean  combinations  of  convex 
polygons  were  used  to  generate  more  complex  structures  (6);  unfortunately.  Boolean 
combinations  can  only  represent  homogeneous  structures.  In  this  report,  we  wish  to  design 
multi-phase  systems,  so  we  must  generalize  the  method  in  reference  6  to  inhomogeneous  systems, 
as  in  reference  8. 


There  are  several  ways  to  generalize  the  Boolean  approach  to  inhomogeneous  structures;  here,  a 
priority-based  approach  is  used.  Simply  stated,  each  terminal  node  in  a  tree  is  assigned  a  priority 
value  (some  real  number,  initially  randomly  generated,  between  1  and  10,  for  example)  and 
material  properties.  Function  nodes  in  this  scheme  are  not  of  any  specific  type,  the  priority 
values  dictate  the  specific  operation;  however,  they  are  always  binary.  To  decode,  a  tree  is 
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traversed  until  a  fully  terminated  subtree  is  reached.  The  priority  values  of  the  two  terminal 
nodes  (at  this  point  they  must  be  convex  polygons)  are  rounded  to  the  nearest  integer  and 
compared,  and  the  polygon  with  the  higher  priority  value  is  overlayed  on  top  of  the  lower 
polygon.  If  the  two  priority  values  are  equal,  then  a  union  operation  is  performed  and  the 
material  values  of  the  node  on  the  left  of  the  tree  are  used.  The  priority  values  of  the  resulting 
polygons  are  preserved  for  future  operations. 

Figure  3  illustrates  the  decoding  process.  Figure  3a  shows  each  point  list  plotted  in  the  plane 
with  their  convex  hulls  as  the  dashed  lines.  The  points  in  each  list  are  represented  as  symbols 
corresponding  to  the  terminal  nodes  in  figure  3b.  Figure  3b  shows  the  tree  structure  defining  the 
priority  combinations.  In  this  example,  material  values  are  represented  with  shading.  Figure  3c 
shows  the  result  of  the  leftmost  subtree.  In  this  case,  the  priority  values  are  equal  when  rounded, 
so  the  result  is  the  union  of  the  two  convex  polygons.  Finally,  figure  3d  shows  the  overlap  of  the 
rightmost  terminal  node  and  the  result  of  the  previous  step.  The  final  result  is  an  example  of  an 
inhomogeneous  structure,  which  cannot  be  encoded  using  only  Boolean  operations.  Now  that  the 
encoding/decoding  scheme  is  defined,  a  set  of  genetic  operators  can  be  defined  to  guide  a 
population  of  such  tree  chromosomes  to  a  solution. 

3.2  Genetic  Operators 

Central  to  any  flavor  of  genetic  algorithm  or  genetic  programming  method,  the  genetic  operators, 
selection,  crossover,  and  mutation  can  be  adapted  to  the  specific  encoding/decoding  scheme.  As 
selection  does  not  involve  chromosome  structure  directly,  standard  techniques  are  used  (16). 
Crossover  involves  hybridizing  two  encoded  chromosomes  by  swapping  subtrees  and  point  lists 
(6)  and  is  discussed  in  detail  in  subsection  3.2.1.  Mutation,  presented  in  subsection  3.2.2, 
increases  genetic  diversity  and  several  different  types  are  implemented.  Many  have  been 
discussed  in  references  6-8,  though  some  have  been  altered  here. 
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(c) 


(d) 


Figure  3.  Example  of  the  decoding  process:  (a)  Point  lists  with  their  convex 

hulls,  (b)  tree  defining  priority  combinations,  (c)  first  subtree  decoded, 
and  (d)  final  result . 
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3.2.1  Crossover 


Crossover  is  the  hybridization  of  two  parent  chromosomes  resulting  in  two  children,  which,  at 
times,  contain  good  traits  from  each  parent.  Typically,  crossover  involves  splicing  parts  of  each 
chromosome  between  parents,  and  possibly  hybridizing  some  nodes.  For  genetic  programming, 
subtrees  are  exchanged  and  the  root  nodes  of  those  subtrees  may  be  hybridized  in  some  way. 
Crossover  for  genetic  programming  can  be  implemented  in  a  straightforward  way;  however,  “tree 
bloat” — the  exponential  growth  of  chromosomes — can  be  a  major  cause  of  population  stagnation 
(17).  A  solution  to  tree  bloat  is  to  design  a  crossover  scheme  that  has  less  variation  than  the 
standard  implementation.  In  other  words,  the  children  should  resemble  their  parents  closely.  In 
reference  8,  this  was  accomplished  two  ways.  First,  a  geometric  (phenotypic)  similarity-based 
crossover  probability  was  used,  and  also  only  terminal  nodes  were  used  in  crossover.  Here,  the 
first  method  is  retained,  while  the  second  is  abandoned.  Only  involving  terminal  nodes  in 
crossover  will  certainly  eliminate  tree  bloat  via  crossover,  though  it  may  be  too  restrictive  for 
more  complicated  topologies.  Instead,  the  population  is  initialized  only  with  single-node  trees 
(i.e.,  a  single  convex  polygon),  and  mutation  is  primarily  responsible  for  altering  the  tree 
structure. 

Selecting  a  chromosome  for  crossover  now  involves  two  steps.  First,  for  a  given  chromosome,  a 
potential  set  of  mates  is  chosen  at  random.  The  number  of  potential  mates  is  set  at 

Nm  =  fiVp/201  ,  (18) 


where  Np  is  the  population  size.  For  each  potential  mate,  a  set  of  random  subtrees  are  chosen 
and  compared.  The  number  of  subtrees  to  compare  for  crossover  is  an  input  parameter  Nca, 
typically  set  to  10.  (Of  course,  in  some  cases  there  may  be  fewer  than  Nc a  unique  combinations 
of  subtrees,  so  only  the  smaller  of  Nca  and  the  total  unique  combinations  are  compared.) 
Regardless,  each  chosen  subtree  is  decoded  and  the  resulting  geometries  are  compared  for 
similarity  using 


(^-Y, 

\  Cl 1  +02/ 


(19) 


where  <i\  and  a2  are  the  areas  of  the  decoded  subtrees,  u\  is  the  area  of  their  intersection,  and  s  is 
some  given  biasing  exponent,  typically  chosen  as  2.5.  The  subtree  combination  with  the  highest 
p0  is  saved  and  used  as  the  similarity  for  that  pair  of  mates.  Once  the  similarity  of  all  mates  is 
computed,  a  crossover  mate  is  chosen  using  roulette  wheel  selection  (16),  and  a  final  crossover 
probability  is  computed 

Pm  =  PcPo,  (20) 


13 


where  pc  is  a  given  crossover  probability,  typically  between  0.9  and  1.  Finally,  a  weighted  coin  is 
flipped  with  probability  pm,  and  if  successful,  the  subtree  pairs  chosen  previously  are  used  for 
crossover.  If  two  terminal  nodes  are  chosen,  their  point  lists  are  hybridized. 

In  addition  to  tree  bloat,  the  point  lists  in  each  terminal  node  are  not  subject  to  a  maximum  size, 
so  they  can  also  grow  exponentially.  To  prevent  this,  point  lists  are  hybridized  using  two-point 
linear  crossover,  and  the  number  of  points  to  swap — while  still  random — is  kept  constant  between 
each  point  list.  Swapping  an  equal  number  of  points  between  two  point  lists  ensures  that  neither 
point  list  changes  size  during  crossover.  Also,  the  points  located  at  each  crossover  site  are 
hybridized  by  taking  a  random,  weighted  sum  of  the  components  of  each  point,  as  is  standard  in 
many  real-coded  genetic  algorithm  implementations.  Priority  values  are  hybridized  in  the  same 
way,  and  currently,  material  values  are  not  changed  as  the  materials  are  chosen  from  a  database. 

3.2.2  Mutation 

The  main  purpose  of  mutation  is  to  insert  new  genetic  information  into  a  population  at  random, 
hopefully  reducing  the  chance  of  premature  convergence.  Given  the  complexity  of  the 
chromosome  used  here  and  its  geometric  meaning,  there  are  many  different  types  of  mutations 
possible.  These  have  been  discussed  extensively  in  references  6-8,  so  only  new  mutations  are 
covered  here.  One  major  difference  in  general  is  that,  wherever  possible,  Gaussian  random 
variables  are  now  used  as  perturbations  rather  than  uniform  random  variables  in  some  specified 
range. 

The  first  new  mutation  is  a  point  addition  mutation,  though  rather  than  randomly  adding  a  point 
anywhere  in  the  valid  region,  the  point  is  added  along  a  line  segment  of  the  convex  hull.  Given 
two  points,  a  point  can  be  generated  lying  at  a  random  point  between  them 

p*  =  thi  +  (1  —  t)h2.  0  <  f  <  1,  (21) 

where  t  is  a  uniform  random  variable  and  h,  and  h,  are  neighboring  points  on  a  convex  hull. 
Adding  p*  to  the  point  list  will  not  affect  the  decoded  shape,  but  it  may  be  useful  in  future 
generations.  In  a  similar  vein,  the  point  addition  mutation  has  been  altered  and  no  longer  simply 
adds  a  random  point.  Point  addition  proceeds  by  duplicating  a  randomly  chosen  point  and 
translating  it  by  a  small  amount.  Previously,  it  was  observed  that  the  point  addition  mutation  was 
destructive,  especially  in  problems  with  a  large  region  size.  Duplicating  and  translating  an 
existing  point  is  less  severe. 

Terminal  node  splitting,  an  important  mutation  operation,  alters  the  topology  of  the  tree  and 
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geometry.  Previously,  terminal  nodes  were  split  in  half  ( 6 )  and  at  the  center  (8).  In  this 
implementation,  a  third  possible  split  is  used,  which  simply  separates  the  terminal  node’s  point 
list  into  two  new  terminal  nodes.  The  split  site  is  chosen  at  random,  and  the  priority  and  material 
values  of  the  original  terminal  node  are  preserved.  This  operation  can  be  useful  for  generating 
new  geometries,  as  the  previous  were  designed  to  generate  an  identical  geometry  when  decoded. 

Finally,  priority  values  and  material  values  can  be  mutated.  Priority  values  are  mutated  by  adding 
a  Gaussian  random  variable  with  a  variance  of  2  and  a  zero  mean.  The  decoded  chromosome 
may  not  change  with  a  change  in  priority  (depending  on  the  priority  values  of  the  other  nodes)  so 
this  can  be  used  with  a  high  rate.  Materials  are  mutated  simply  by  randomly  choosing  a  new 
material  from  the  database. 

In  its  entirety,  the  genetic  programming  method  presented  here  is  similar  to  references  6-8,  but 
with  a  few  important  changes.  These  changes  lead  to  a  more  efficient  and  effective  method  as  is 
demonstrated  in  section  4. 


4.  Results 


The  following  are  a  few  examples  of  the  method  described  previously  using  different  materials. 
The  following  examples  were  run  on  a  parallel  computing  system,  each  using  20,  four-core 
computers  giving  a  total  of  80  nodes.  Input  parameters  for  each  example  are  given  in  table  1. 

4.1  Example  1 

In  the  first  example,  two  different  materials  were  available  to  the  algorithm:  epoxy  and  lead,  with 
material  properties  taken  from  reference  2  and  epoxy  fixed  as  the  matrix.  A  square  unit  cell  was 
used,  with  side  lengths  of  1  m.  The  best  result  found,  shown  in  figure  4  with  light  gray  indicating 
epoxy  and  dark  gray  indicating  lead,  has  a  relative  bandgap  size  of  0.68  and  required  260 
generations.  Figure  5  gives  the  band  diagram  for  this  example,  and  the  center  frequency  of  the 
bandgap  is  about  850  Hz. 
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Table  1.  Genetic  programming  input 
parameters. 


Parameter 

Value 

Size 

300 

Tournament  Size 

4 

Population 

Crossover  Rate 

85% 

Elitism 

on 

Delete  nodes 

0.5% 

Delete  one  point 

2% 

Add  duplicate  point 

0.5% 

Split  hull  segment 

0.5% 

Geometry  transform 

1% 

Priority 

1% 

Material 

1% 

Mutation 

Prune 

5% 

Hull 

1% 

Point  translation 

1% 

Split  in  half 

1% 

Split  a  hole 

1% 

Split  list 

1% 

Combine 

6% 
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Figure  4.  Results  of  example  1;  light  gray  indicates  epoxy  and  dark  gray  indicates 
lead. 


Figure  5.  Band  diagram  for  example  1. 
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4.2  Example  2 


Next,  a  third  material  was  added,  silicone  rubber,  again  with  material  properties  taken  from 
reference  2;  however,  the  algorithm  was  allowed  to  choose  the  background  material.  These 
materials  were  chosen  in  reference  2  to  construct  a  low  frequency,  local  resonance  based  bandgap, 
with  a  design  of  a  silicone  rubber-coated  lead  sphere  embedded  in  an  epoxy  matrix  laid  out  in  a 
cubic  arrangement.  Again,  a  square  unit  cell  was  used  with  side  lengths  of  1  m.  A  bandgap  size 
of  1.56  was  achieved  after  565  generations.  The  best  result  found  is  shown  in  figure  6.  In 
contrast  to  the  design  in  reference  2,  the  genetic  programming-generated  design  does  not  coat  the 
entire  mass  of  lead  in  silicone  rubber,  though  it  still  shares  common  features  such  as  an  epoxy 
background  with  lead  embedded  in  silicone  rubber.  The  band  diagram  for  this  design  is  shown  in 
figure  7,  which  shows  the  center  frequency  of  the  band  is  around  550  Hz. 
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Figure  6.  Results  of  example  2;  light  gray  indicates  epoxy,  dark  gray  indicates 
lead,  and  black  represents  silicone  rubber. 
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Figure  7.  Band  diagram  for  example  2. 
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5.  Conclusions 


A  genetic  programming  method  was  applied  to  design  the  inclusion  of  phononic  bandgap 
materials  to  maximize  the  bandgap  size.  The  geometry  was  encoded  using  a  tree  structure  that 
combines  convex  polygons  using  a  priority-based  overlapping  scheme.  The  algorithm  could 
choose  from  a  database  of  materials  for  both  the  inclusion  and  the  matrix.  Results  were  given  for 
three  materials:  lead,  silicone  rubber,  and  epoxy.  One  example  restricted  the  method  to  use  only 
lead  and  epoxy  with  epoxy  fixed  as  the  matrix  material.  In  another  example,  all  three  materials 
were  allowed  with  no  fixed  matrix.  Bandgaps  were  achieved  for  each  material  system. 

In  the  future,  the  method  could  be  extended  to  a  multi-objective  problem,  where  bandgap  size 
could  be  simultaneously  optimized  with  another  goal  such  as  center  frequency  or  inclusion  size. 
Additional  materials  could  be  used  for  specific  applications. 
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